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Field of the Invention 

[0003] Generally speaking, the present invention relates to systems and methods that model 
20 the impact of a medium on the appearances of encompassed light sources. More particularly, 
the present invention relates to systems and methods which model multiple scattering of light 
from a light source encompassed by the medium, and that may then determine characteristics 
of the medium or alter images of the light source to remove or add effects of the medium. 

25 Background of the Invention 

[0004] Arguably, nothing in the world impacts the daily lives of people more than the 
weather. 

Perhaps nowhere is the effect of weather more strongly felt than on the runways of world's 
airports. Haze, mist, fog, and rain can and do prevent airplanes from being able to take off 
30 and land. Not only do these weather conditions present an inconvenience to travelers in the 
form travel delays, these conditions can be life threatening when fog, for example, makes it 
impossible for an airline pilot to see the runway. Naturally, such weather conditions also 
impact life in countless other ways as well. 

[0005] Because of these effects, researchers have attempted to build machines and develop 
35 models that can accurately gauge the weather. For example, at many airports, visibility 
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meters have been put in place for determining the visibility at discrete locations around the 
grounds of the airport. Such machines detect the impact of atmospheric conditions (such as 
fog) on a beam of light transmitted from a light source, such as a laser, to a detector. These 
machines rely on the principle that particles in the atmosphere scatter light incident on those 
5 particles. Such machines are problematic, however, in that they only sample light in a small 
region (e.g., a two foot by two foot region) and are based on single scattering. Moreover, 
because of the complexity of these machines, they are quite expensive. 
[0006] In contrast to the single scattering principle used in these machines is the principle 
of multiple scattering. One easily understood illustration of multiple scattering is the manner 

10 in which bad weather impacts the appearance of light sources, especially when viewed at 
night. This impact is typically seen in the form of a glow or halo that appears to surround a 
light source. The appearance of this glow or halo is the result of light rays being deflected 
multiple times by particles in the atmosphere as they leave the light source, such that the rays 
appear to be originating from an area surrounding the light source. 

15 [0007] Attempts have been made to model the effects of multiple scattering on light rays in 
the atmosphere. For example, one approach has been to assume that the paths of light 
traversal are random and then to apply numerical Monte-Carlo techniques for ray tracing. 
Monte-Carlo techniques are problematic, however, in that they attempt to predict the paths of 
light based upon each photon impacting each particle in the medium. Such an approach is 

20 highly dependent upon the specific properties of the medium, such as the density of the 

medium, type of particles in the medium, the size of the medium, etc. As a consequence of 
this, millions of rays must be traced through the atmosphere to accurately model the multiple 
scattering using Monte-Carlo techniques. Thus, when attempting to model a single light 
source in a single image, this technique can take several hours using current computer 

25 systems. Obviously, such a long processing time is inadequate for real-time or near-real-time 
image processing. 

[0008] Monte-Carlo techniques are also problematic in that they do not converge for pure 
scattering media. This may be unacceptable when trying to model multiple scattering of 
certain media. 

30 [0009] Another technique for modeling the impact of particles in the atmosphere on light is 

to use the plane parallel model. This model has been used in fields such as atmospheric 

optics and astronomy in which the medium being observed (e.g., the atmosphere around the 

earth) is illuminated by a distant light source such as the sun or moon. The plane parallel 

technique is problematic, however, in that it assumes that the light source is entering the 
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medium from an infinitely far away distance, which is not the case when modeling the impact 
of an atmosphere on the appearance of a light source that is inside that atmosphere. 
[0010] Accordingly, other techniques for modeling the multiple scattering of light from 
encompassed light sources and for using those models are desired. 

5 

Summary of the Invention 

[0011] The present invention provides systems and methods for modeling the impact of a 
medium on the appearance of an encompassed light sources using a Legendre polynomial 
series solution to a Radiative Transfer Equation for Spherical Media (RTE-SM) called an 

10 Atmospheric Point Spread Function (APSF). Using this APSF, it is possible to determine 
characteristics of the medium causing the multiple scattering of the light from the 
encompassed light source. For example, by observing a street light in bad weather at night, 
using the APSF, it is possible to determine whether the bad weather is haze, mist, fog, or rain. 
Similarly, the APSF may be used to estimate the size of particles in a liquid. It is also 

15 possible using the APSF to remove and/or add an effect of the medium on a light source 
captured in an image. 

[0012] In accordance with an embodiment of the invention, a system features an image 
acquisition device that captures an image of a light source encompassed in a medium and a 
computer that receives the image from the image acquisition device, identifies the light 
20 source in the image, models the multiple scattering of light from the light source in the 
medium using a RTE-SM, and determines one or more properties of the medium using the 
RTE-SM. 

[0013] In accordance with another embodiment of the invention, a system features a 
computer that receives an image, identifies a light source in the image, generates a model of 
25 multiple scattering of light from the light source in a medium in the image using a RTE-SM, 
and alters the image based upon the model. 

[0014] In accordance with another embodiment of the invention, a method features 
capturing an image of a light source, identifying the light source in the image, modeling the 
multiple scattering of light from the light source in the medium using a RTE-SM, and 
30 determining one or more properties of the medium using the RTE-SM. 

[0015] In accordance with another embodiment of the invention, a method features 
receiving an image, identifying a light source in the image, generating a model of multiple 
scattering of light from the light source in a medium using a RTE-SM, and altering the image 
based upon the model. 
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[0016] In accordance with another embodiment of the invention, a method of monitoring 
weather conditions features locating an image acquisition device in a location suitable for 
capturing images of multiple light sources, sequentially aiming the direction of the image 
acquisition device at each of the multiple light sources, capturing an image of each of the 
5 multiple light sources, identifying the light source in each of the images, modeling multiple 
scattering of light from the light source using a RTE-SM, and determining at least one of the 
forwarding scattering parameter, the optical thickness, the visibility of the atmosphere around 
the light sources using the RTE-SM. 

[0017] In accordance with another embodiment of the invention, a method of monitoring 
10 weather conditions features locating a first image acquisition device in a first location 
suitable for capturing images of a first light source, locating a second image acquisition 
device in a second location suitable for capturing images of a second light source, capturing 
an image of each of the first light source and the second light source, identifying the light 
source in each of the images, modeling multiple scattering of light from the each light source 
15 using a Radiative Transfer Equation for Spherical Media, and determining at least one of the 
forward scattering parameter, the optical thickness, the visibility of the area using the 
Radiative Transfer Equation for Spherical Media. 

[0018] In accordance with another embodiment of the invention, a method of monitoring 
weather conditions features locating an image acquisition device in a location suitable for 

20 capturing images of a light source, capturing multiple images of the light source, averaging 
the images to produce an averaged image, identifying the light source in the averaged image, 
modeling multiple scattering of light from the light source as captured in the averaged image 
using a Radiative Transfer Equation for Spherical Media, and determining at least one of the 
forward scattering parameter, the optical thickness, the visibility of an area using the 

25 Radiative Transfer Equation for Spherical Media. 

Brief Description of the Drawings 

[0019] FIG. 1 shows four images illustrating multiple scattering of light from a light source. 
[0020] FIG. 2 shows an illustration of a prior art approach to modeling light scattering. 
30 [0021] FIG. 3 shows an illustration of an embodiment of an approach to spherical modeling 
in accordance with the present invention. 

[0022] FIG. 4 shows another illustration of an embodiment of an approach to spherical 
modeling in accordance with the present invention. 
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[0023] FIG. 5 shows an illustration of light incident on and exiting from an infinitesimal 
element of the atmosphere in accordance with an embodiment of the present invention. 
[0024] FIG. 6 shows an illustration of classifications of different ranges of forward 
scattering parameters q in accordance with an embodiment of the present invention. 
5 [0025] FIG. 7 shows a graph of an Atmospheric Point Spread Function (APSF) in 
accordance with one embodiment of the present invention. 

[0026] FIG. 8 shows another graph of an APSF in accordance with one embodiment of the 
present invention. 

[0027] FIG. 9 shows a graph of the coefficient term m value versus optical thickness T in 
10 accordance with one embodiment of the present invention. 

[0028] FIG. 10 shows a graph of an APSF in which ringing is present because not enough 
coefficient terms m were used in calculating the APSF in accordance with one embodiment 
of the present invention. 

[0029] FIG. 1 1 shows a flowchart illustrating a process for calculating the forward 
15 scattering parameter q and the optical thickness T in accordance with one embodiment of the 
present invention. 

[0030] FIG. 12 shows a flowchart illustrating a process for altering an image in accordance 
with one embodiment of the present invention. 

[0031] FIG. 13 shows an illustration of two possible arrangements of cameras and mediums 
20 being imaged in accordance with an embodiment of the present invention. 

[0032] FIG. 14 shows an illustration of another possible arrangement of a camera and a 

medium being imaged in accordance with an embodiment of the present invention. 

[0033] FIG. 15 shows an illustration of an atypical shape and size light source in 

accordance with one embodiment of the present invention. 
25 [0034] FIG. 16 shows a block diagram of equipment in one embodiment of the present 

invention. 

[0035] FIG. 17 shows a block diagram of equipment in another embodiment of the present 
invention. 

30 Description of Particular Embodiments of the Invention 

[0036] Turning first to FIG. 1, an example of a street lamp in four different atmospheric 

conditions is illustrated in images 25, 26, 27, and 28. As can be seen, in the first image 25, 

no glow surrounds the individual globes of the street lamp. This is the case because the 

weather as shown in the image is clear. In image 26, however, it is apparent that the street 
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lamp is now encompassed by mild fog and a slight glow appears around the individual globes 
of the lamp. Next, in image 27, the individual globes of the street lamp are all but obscured 
by the dense mist that now encompasses the street lamp. Finally, in image 28, the individual 
globes of the street lamp cannot be discerned because of the highly dense haze encompassing 
5 the street lamp. The glows surrounding the lamp in images 26, 27, and 28 are the result of 
multiple scattering of the light rays that are exiting the individual globes of the lamp. 
[0037] Referring to FIGS. 2 and 3, two different scenarios for light from a light source 
impacting a medium are illustrated. First, as shown in FIG. 2, a light source 20 impacts a 
homogenous medium 21 from an infinite distance. Light rays 22 from source 20 travel 

10 through medium 21 and scatter at point 23 in various directions as indicated by scattered light 
rays 24. Because the light rays from source 20 impact the plane of medium 21 from an 
infinite distance in a parallel fashion, this is referred to as a plane parallel medium. Such a 
medium is not a good representation of an isotropic light source within an encompassing 
medium (like the street lamp in FIG. 1) because, unlike the parallel impact of the rays from 

15 the external light source in the plane parallel medium, light rays from a light source within an 
encompassing medium enter the medium radially from the light source. 
[0038] A better representation is the spherical medium illustrated in FIG. 3. As shown, a 
light source 30 is encompassed by a homogenous medium 31. Light source 30 might be a 
street lamp within a medium 31 such as dense fog. Light rays 32 from source 30 radiate 

20 isotropically from source 30. Some of the rays 32 reach a "pin hole camera" 33 at a radius 35 
from source 30 and then impact a surface within the camera at different points 34 and at 
different angles 0, 36, relative to and normal to the camera's pin hole. In the spherical 
medium representation, it is assumed that the medium is homogeneous and infinite in extent. 
[0039] An alternate view of the spherical medium 31 shown in FIG. 3 is illustrated in FIG. 

25 4. As can be seen, a medium 41 encompasses a light source 40 that is outputting a light ray 
43. Light ray 43 scatters at points (particles) 47 and 48 and passes through pin hole 44 and 
onto camera plane 45. The angle between the last leg of ray 43 hitting plane 45 and the 
straight line 42 between light source 40 and plane 45 is represented by angle 6, 46. Other rays 
that reach plane 45 from source 40 through pin hole 44 may hit plane 45 at other angles also 

30 represented by 
angle 6. 

[0040] Generally speaking, light scattering in a spherical medium, such as the earth's 
atmosphere, water solutions, blood, tissue, smoke, milk, etc., can by modeled by the physics- 
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based theory of radiative transfer. As illustrated in FIG. 5, the key idea in this approach is to 
investigate the difference between light 50 incident on and light 52 exiting from an 
infinitesimal volume 51 of the medium. The change in flux through such volume 51 can then 
be given by an integro-differentiation equation called a Radiative Transfer Equation for 
5 Spherical Media (RTE-SM). 

[0041] More specifically, a solution to the RTE-SM can be determined as follows. First, 
the directional distribution of light incident on a particle, called the scattering or phase 
function of the particle (P(cosa)), must be considered. For most atmospheric conditions, this 
phase function is known to be symmetric about the direction of incident light. Thus, the 
10 phase function depends only on the angle a between the directions of incident and scattered 
radiation. The cosine of this angle is determined as: 

cosa = V(l -\i 2 )(\-\i a ) costy -<?)) (1) 
wherein: 6,(() is the direction of the scattered radiance; 

8\<t>' is the direction of the incident irradiance; 
15 |X = cos9;and 

\l* = cos0\ 

[0042] Axially symmetric phase functions (both isotropic and anisotropic) can be expanded 
using Legendre polynomials L k as: 

P(cosar) = ]T W k L k (cosa) (2) 

k=0 

20 [0043] The exact shape of the phase function depends on the size of the scattering particles. 
For instance phase functions of small particles (e.g., air molecules) have a small peak in the 
direction of incidence. Isotropic and Rayleigh phase functions describe the scattering for air 
molecules. An Isotropic phase function that may be used in the present invention can be 
expressed as: 

25 P(cosa) = W 0 (3) 

wherein: W k =0 for all k>=l . 

A Rayleigh phase function that may be used in the present invention can be expressed as: 

3W 

P(cosa) = — 2.(1 + cos 2 a) (4) 
4 

wherein: Wi=0; 
30 W 2 =Wo/2; and 

W k =0 for all k>=3. 
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Yet another phase function that may be used in the present invention can be expressed as: 
P(cosa) = W 0 (l + xcosa) (5) 

wherein: -1 <= x <=1; 

W,=Wox; and 
5 W k =0 for all k>=2. 

Still another phase function that holds for particles of various sizes, that may be used in the 
present invention, and that is used herein for purposes of illustration is the Henyey-Greenstein 
phase function: 

P(cos a) = l-q 2 /(l-Kj 2 -2qcosa) 3/2 (6) 
10 wherein: q is the forward scattering parameter and is between 0 (i.e., isotropic 
scattering) and 1 (i.e., anisotropic scattering - all scattering in forward 
direction); and 

the coefficients of the Legendre polynomial expansion are: 
W k =W 0 (2k+l)q k for k>=l. 
15 [0044] Note that the W 0 coefficient is called the single scattering albedo, which is the 

fraction of flux scattered. When W 0 =l, the medium is purely scattering, and if W 0 <1, there is 
absorption. 

[0045] Next, the Radiative Transfer Equation for Spherical Media (RTE-SM) for a 
spherically symmetric atmosphere may be written as: 
-w 1 _ // 2 i 2ff+l 

20 M J£. + L-£- = - IiTyM)+ jjp(cosa)I(r,M'W^' O) 

oT T 4tt 0 _j 

[0046] Finally, by expanding the Henyey-Greenstein phase function in terms of Legendre 
polynomials, these equations can be solved as the Atmospheric Point Spread Function 
(APSF) (the derivation for this equation is presented at the end of this Detailed Description): 

I(T, = £ (g m (T) + gm + WUbi) (8) 

m=0 

25 wherein: I(T,|ll) is the intensity of light; 

Lm(H) is the Legendre polynomial of order m, 

i a m . 2 



a,n=m+l; 
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q is the forward scattering parameter; 
T is the optical thickness of the medium, 

V 

5 o is the extinction coefficient which denotes the fraction of flux lost due to 

scattering within a unit volume of the medium; 
V is the visibility of the medium; 

R is the distance from the light source to the camera's pin hole; and 
|l = cos8. 

10 [0047] In these equations, the function g m (T) captures the attenuation of light in the 
medium, whereas the Legendre polynomial Lm([i) explains the angular spread of the 
brightness observed due to multiple scattering. This model is valid for both isotropic (i.e., 
q=0) as well as anisotropic (i.e., 0<q<=l) scattering, and thus describes glows under several 
weather conditions. The series solution in equation (8) does not converge for T<=1. 

15 Nevertheless, multiple scattering is very small for small optical thicknesses (i.e., when T<=1) 
and thus glows are not seen. 

[0048] As illustrated in FIG. 6, in an embodiment in which the APSF is being used to 
model multiple scattering in the earth's atmosphere, different values of q, 60, can be used to 
represent different weather conditions 61. For example, values of q, 60, ranging from 0.0 to 

20 0.2 may represent "air," values ranging from 0.2 to 0.7 may represent "small aerosols," 
values ranging from 0.7 to 0.8 may represent "haze," values ranging from 0.8 to 0.9 may 
represent "mist" or "fog," and values ranging from 0.9 to 1.0 may represent "rain." Thus, as 
can be seen, q is an indicator of the size of the particles causing the multiple scattering. It is 
apparent that the boundary values in this example, e.g., 0.7, are approximate and do not 

25 strongly indicate which of the two bordering weather conditions is present. 

[0049] Turning to FIGS. 7 and 8, illustrations of the results of the APSF are shown for two 
different extinction coefficients, o. These results are normalized so that the value of intensity 
I (on the y-axes 70 and 80) ranges from 0 to 1. On the x-axes 74 and 84, the values 
represented are the angles of 9, 46 (FIG. 4), from the camera's pin hole 44 (FIG. 4) relative to 

30 the camera's plane 45 (FIG. 4) for the corresponding values of I. As can be seen in FIG. 7, 
the curves 71, 72, and 73 represent the APSF for the spherical model for medium when 
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containing "small aerosols," "haze," and "fog," respectively. As also shown, the values for q 
for these three curves are approximately 0.2, 0.75, and 0.9, respectively. Turning to FIG. 8, it 
can be seen that curves 81, 82, and 83 similarly represent the APSF for the spherical model 
for a medium when containing "aerosols," "haze," and "fog," respectively. Also similarly, 
5 the values for q for these three curves are approximately 0.2, 0.75, and 0.9, respectively. 
[0050] As can be seen from these curves, larger particles (such as fog particles) result in 
greater forward scattering, and thus will result in narrower or more pointed APSF curves. As 
can also be seen from the curves, curves 81, 82, and 83 in FIG. 8 are wider than the 
corresponding curves 71, 72, and 73 in FIG. 7. This is because the extinction coefficient (a) 

10 (and hence the optical thickness (T)) is much greater in FIG. 8 than in FIG. 7 as a result of the 
atmosphere in FIG. 8 being much denser than that in FIG. 7. 
[0051] Before equation (8) can be solved, it is necessary to determine how many 
coefficients m to use. As can be seen in FIG. 9, the number of coefficient terms m, 90, to use 
is a function of the optical thickness T, 91, of the medium. Because T may not be known 

15 when solving the APSF, a value of m greater than 200 should be sufficient to solve the APSF. 
In the event that too few coefficient terms m are used, ringing will appear on the APSF curve 
such as that shown in curve 100 of FIG. 10 at points 102. Note that in FIG. 10, when m=200 
in curve 101, no ringing is apparent. 

[0052] Using equation (8), it is now possible to determine the weather conditions and 
20 visibility from an image or to alter the glow (if any) in an image based upon desired weather 
conditions and visibility. In order to determine the weather conditions, the process of FIG. 1 1 
may be used. As shown, at step 102, an image may be captured using any suitable image 
acquisition device, for example a camera. Next, at step 103, the image is transferred to a 
suitable computer. The process then determines the values of Iq (the intensity of the source), 
25 R (the range from the source to the camera), and m (the number of coefficient terms to use) at 
step 104. In order to determine Io when I 0 is captured by more than one pixel in the CCD of 
the camera, it may be necessary to treat Io as a light source with an atypical shape or size, as 
discussed below in connection with FIG. 15, and use a thresholding technique to identify the 
shape function S of the source and then use the shape function to convolve the intensity 
30 values Io. At step 105, the computer may then average the intensity values I captured by the 
camera along the radial contours of the image. In order to increase the signal-to-noise ratio 
of the image (e.g., when using an inexpensive camera as the image acquisition device), in 
some embodiments the average intensity values of intensity I and the value of source 
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intensity Io may also be individually averaged over time. Next, using a suitable curve fitting 
algorithm, such as multidimensional unconstrained nonlinear minimization (i.e., Nelder- 
Mead search), the values of T and q can be calculated so that equation (8) approximates the 
curve provided by the measured values of I at step 106. These values of T and q can then be 
5 plugged into equation (8) to determine if any ringing exists in the resultant APSF curve as 
shown in FIG. 10. In the event that any ringing is present, step 108 may then be performed to 
increase the value of m, after which step 106 is repeated. Otherwise, if no ringing is present, 
the weather in the image may be looked-up from the determined value of q at step 109. The 
visibility may then be determined at step 1 10 based upon the values of T and R. Finally, at 

10 step 1 1 1, the data is output to the desired user or application. 

[0053] Although the process in FIG. 1 1 used a known value of R to calculate values of T 
and q, it is possible to use known values of T to calculate values of R and q, and to use 
known values of q to calculate values of T and R in substantially the same manner. 
[0054] In order to alter glow (if any) in an image based upon desired weather conditions 

15 and visibility, the process in FIG. 12 may be performed. First, at step 120, the image to be 
altered may be selected. Next, at step 121, the light source in the image may be determined. 
The values of Io (the intensity of the light source in the image), R (the range to the light 
source in the image), and m (the number of coefficient terms to use) may be determined at 
step 122. In order to determine Io when Io is captured by more than one pixel in the CCD of 

20 the camera, it may be necessary to treat I 0 as a light source with an atypical shape or size, as 
discussed below in connection with FIG. 15, and use a thresholding technique to identify the 
shape function S of the source and then use the shape function to convolve the intensity 
values Io. Next, at step 123, the desired values of T and q are determined. Using equation 
(8), the APSF is then determined at step 124. This determination provides the values of I for 

25 the different values of 8. At step 125, the process then determines if there is any ringing in 
the values of I as shown in FIG. 10. If ringing is present, the number of coefficient terms m 
is then increased at step 126 and step 124 repeated. Otherwise, if no ringing is present, the 
image is then convolved or deconvolved using the APSF at step 127 to add or remove a glow. 
Finally, at step 128, the image is output to the desired user or application. 

30 [0055] In the event that the lens optics of the camera used in embodiments of the present 
invention has a small field of view (e.g., when using a zoom lens), intensities (I) from some 
directions may be cut-off. In such case, the glow surrounding the light source can be 
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approximated by multiplying the intensity values at the outer edge of the glow by a constant 
attenuation factor (e" T /R 2 ) and then using the product of that multiplication for all values of I. 
[0056] Alternate equations for the definition of \i may be necessary to take into account the 
geometry of any particular arrangement of camera and medium being imaged. For example, 
5 as illustrated in FIG. 13, when an orthographic camera 136 is used instead of a pin hole 
camera 138, angle 9 in the equation \i=cos 9 should be equal to angle 132 relative to the 
perpendicular of plane 130 rather than angle 134. As another example, as illustrated in FIG. 
14, when a pin hole camera is used that is not within the medium such as in arrangement 140, 
the variable |i may be represented by the following: 



[0057] In certain instances, it may be desirable to model a light source of an atypical shape 
or size, such as the neon sign 150 illustrated in FIG. 15. Such sources can be assumed to be 
made up of several isotropic source elements with varying radiant intensities Io(x,y). By the 
principle of superposition of light (it is assumed that the source elements are incoherent), the 
15 intensities due to different source elements can be added to yield the total intensity due to the 
light source. It is also assumed that the entire area of the light source is at the same depth 
from the observer and passes through the same atmosphere. The image of a light source of 
arbitrary shape and size may then be written as the following convolution: 



Because APSF is rotationally symmetric, this two-dimensional convolution can be replaced 

by two one-dimensional convolutions in order to render the source more quickly. 

[0058] In order to determine the shape function S, it may be necessary in certain cases to 

25 determine the shape function by performing thresholding or high-pass filtering on the 

corresponding area of the image. Another approach is to determine the APSF from a nearby 
first light source and then use that APSF to deconvolve the image of the second source and 
remove the glow around the second source. This may be advantageous where the shape of 
the first light source is simple, or it is likely to be easier to perform thresholding on that light 

30 source. For example, taking a first source and applying a simple thresholding function to the 
source may give a shape function Si for that source. Assuming that the radiant intensity of 



10 




20 wherein: 



I=(IoS)*APSF 

S is a characteristic shape function that is constant over the extent of the light 
source (not including the glow). 
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the first source is constant across the source, the APSF can then be recovered from the 
following function: 

APSF=I*(IoSi)° 

This normalized APSF can then be used to deconvolve the image of the second source and 
5 remove the glow around the second source using the following function: 
S 2 =I*(I 0 APSF)- 1 

[0059] With color images, it is necessary to calculate q and T separately for each of the 
three colors making up the images. In doing so, each value of I and I 0 should correspond to 
the appropriate color. 

10 [0060] Turning now to FIG. 16, one embodiment of equipment for implementing the 
present invention is shown. As can be seen, an image acquisition device 162 is used to 
capture an image of a lamp 160 possibly having a glow 161. The field of view of device 162 
is represented by a field 165. The distance from device 162 to lamp 160 is represented by a 
line 166. Once an image is captured, the image may then be transferred to a computer 163 

15 using a link 164. Device 162 may be any suitable image acquisition device, such as a Kodak 
DCS 760 camera or any other suitable camera, video camera, optical sensor, etc., that is 
radiometrically calibrated. The computer may be any suitable computer, such as a personal 
computer, a micro computer, a server, a workstation, a mainframe, a laptop, a PDA, a 
computer incorporated in a telephone (e.g., a cellular telephone), a computer in device 162 or 

20 another camera, etc. Link 164 may be any suitable link such as a USB interface (e.g., 
versions 1.0 or 2.0), IEEE 1394 interface, a serial link, or any other suitable link. An 
alternate embodiment is shown in FIG. 17. As illustrated, camera 170 may be a film or 
instant photography camera. An image from that camera may then be placed in any suitable 
scanner 172 for transferring the image into computer 163. In such an implementation, the 

25 aperture settings of camera 170 may need to be fixed to match a radiometric calibration of 
scanner 172. 

[0061] In both of these embodiments it may be necessary to take into account the Optical 
Transfer Function (OTF) of any camera or scanner used in accordance with the invention. 
First, it may be desirable to check the camera or scanner to verify that the OTF of the camera 
30 or scanner does not affect the frequencies of the APSF. This can be accomplished by (i) 
plotting all possible frequencies of the APSF by taking the Fourier transform of I(T,ji) in 
equation (8) for all possible values of T (e.g., by selecting maximum R and V values and 
using the formula for T set forth above) and \i, and (ii) measuring the OTF using resolution 
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charts. If the OTF is high (e.g., close to 1) for the main frequencies of the APSF, then the 
lens is acceptable. Otherwise, it may be desirable to select a better lens. Second, it may be 
desirable to deconvolve the images taken with the camera or scanner using the camera's or 
scanner's OTF. 

[0062] In accordance with the present invention, the equipment illustrated in FIGS. 16 and 
17 may be configured at any suitable location and arranged so that a single image acquisition 
device that is controlled by a motor 167 is used to monitor many different locations. For 
example, at an airport, the image acquisition device may be located on the control tower so 
that it can be aimed at various lights throughout the grounds of the airport. The computer 
may then be configured to control the motor so that the image acquisition device is 
periodically directed at each of these lights, an image captured, and the processing in FIG. 1 1 
performed to determine the values of T and q in each of the various directions. If the 
appropriate image acquisition device is selected, the lights may be miles away from the 
camera. Such an arrangement could also be used on highways or in any other suitable 
location. 

[0063] Similarly, multiple image acquisition devices could be used and each could be 
focused on a unique one or set of light sources. The images captured by those devices could 
then be processed as described above. 

Derivation of Equation (8) 

[0064] In this section, the derivation of an analytic form for I(T,|i) is provided. First, 
generally speaking, the partial derivative I with respect to |i is eliminated by integrating the 
RTE-SM with respect to \i over the range [-1,+1]. Here, the partial derivative is eliminated 
by using exact integrals instead of approximating the integrals by linear summations. 
[0065] More particularly, by integrating the phase function P over the azimuth angle <|>\ a 
function that does not depend on the azimuth angle can be defined as: 



and, n=cos 0 and n'=cos 9'. As we shall see, the use of P (0) simplifies the mathematics 
involved in modeling the multiple scattering around a point light source. Substituting into the 
following equation for an RTE-SM: 




(9) 



where, 



cos a = ////'+V(1 - M 2 K 1 " M' 2 ) cos(0 - § ) 



(10) 
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1 — a 2 31 1 2 * +1 

^^+-^-^ = -1^,//)+— J|PC" > <z>;// , ^')l(T,// , )9// , af (11) 

we get: 

//^ + i^^ = -I(T,//) + i|p (0 >Ct/,^)I(T,// , )a//'. (12) 
A function Q m (|J.) for some m>0 is defined such that: 

!. w= - M^M>) oQ . M= i^L (13) 

3// m(m + 1) 



Consider the integral: 



fd-)" 2 )Q n 



-1 



dM (14) 



Integrating by parts and using (13), it has been shown that the partial derivative with respect 
to \i can be eliminated: 



91 



10 jd-// 2 )Qj y- d// = fl(T,//)L m 8// . (15) 



When there is no confusion, the parameters \i and T can be dropped for brevity. Multiplying 
(12) by Q m and integrating with respect to |X over [-1,+1], provides: 

|//Q m ^d M + J^Q m |^ = - |Q m i^ + \ \qJm |p (0) (M^WT^W . (16) 
Substituting equation (15), the RTE can be written as: 

15 J/<Q n + ^ jL m ia// = - |Q m ia// + i jQ m (//)3//JP (0) (//,// , )I(T,// , )a^'. (17) 
-1 " 1 A -1 -1 ^ -1 -1 

Thus the partial derivative with respect to (i is eliminated from the RTE-SM leaving only the 
derivative with respect to the optical thickness T. 

[0066] Next, it is assumed that a solution I m (T,^) to (17) is a product of two functions - 
g(T) depending only on the optical thickness T, and f(|i) expressing the angular dependence. 
20 Mathematically, this is represented as: 

I m (T,//) = g m (T)f m (//) (18) 
Substituting into (17), provides: 

g' m )mQJJm+^\lJJM + g m + \QJjM + fQ m d/fP (0) WytJ/tW = 0 (19) 
-1 1 -1 -1 1 -I -1 
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Since the RTE-SM is spherically symmetric (does not depend on the azimuth angle <()), a 
Legendre polynomial expansion (and not a spherical harmonic expansion) of I is possible. 
Therefore, suppose 

U/O^^+L. (20) 
5 for some m>0. As shall be shown, the above form is key to get a non-trivial solution to the 
RTE-SM ensuring none of the terms of the RTE-SM go to zero. Since the phase function P is 
symmetric about the direction of incident light, P can be expanded using Legendre 
polynomials as 



P(cosa) = ]T W k L k (cosa) (21) 

10 Then, by using (10), and the following identity: 

L k (cos 0) = L k (//)L k (//') + 2]T L n k (//)L n k (//')cos n(<f> - fi) (22) 

the following is provided: 

P (0) (//,//') = £ W k L k (//)L k (//') (23) 

Similarly L\(\i) and |iL' k (^i) can be expanded using a Legendre polynomial series: 
15 L' k (//) = (2k - l)L k _, (//) + (2k - 5)L k _ 3 (//) + ... 

{JL\ (//) = kL k (//) + (2k - 3)L k _ 2 (//) + ... (24) 
Substituting equations (13), (20), (23), and (24) into equation (19) and simplifying each term 
using the orthogonality of Legendre polynomials: 

2 



[ULdju = if i=j=n; 

J 1 J ^ 2n + l 



20 jhiL-dju = 0 otherwise; 

-i 

[0067] The above property can be used to greatly simplify the mathematics as detailed 
below. 

Term 1 in equation (19): 

+i , +i 

-= SLiS f///T' \rt -i-T 



_ J , m(m + l). J , 

25 =— ^-|{mL m +(2m-3)L ro _ 2 + ...}(L m _ 1 +L m 

m(m + 1) _j 
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+1 



= d!f^J {mLJ(Lja/ ' 

m(m + 1) _j 



=g', 



r 2 i 



Term 2 in equation (19): 



_g 



Term 3 in equation (19): 



4 m(m + l)j; 



(25) 



(26) 



10 



= _ TT7: j((2m-l)L m _, +...}(L m . 1 +L m )d/l 
m(m + 1) _ J , 



2g„ 



m(m + 1) 



Term 4 in equation (19): 



+1 +1 



-1 -1 



(27) 



gr 



2m(m 



— \h' m {n)dM |j] W k L k (//)L k in') j(L m _, (//') + L m (/Otf// 



2m(m • 



15 



2m(m- 



= is f^±l(2m - 1) 

2m(m + l)V2m-l ) 2m -1 



2g 



m(m 



+ l)Um-lJ 



(28) 



[0068] The fact that P(0) can be expressed as products of Legendre polynomials is clearly a 
considerable advantage in simplifying the terms of equation (19). Substituting the 4 terms 
20 from equations (25), (26), (27), and (28) into equation (19), provides: 
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g 



I §m f 2 ^ ^ 2g m 2g m ( W mM V (29) 

+ 1)J T Um + lJ m(m + l) m(m + l)Um-lJ 



^(m + l)(2m- 

The above equation can be simplified and written in a concise manner as: 



g'.+Y«-+g-A=0 (30) 



wherein, 



a m - m + 1 

^ m ^ m J{ 2m — 1, 
For the Henyey-Greenstein phase function: 



r2m+lY,_ W.. 



The solution to (30) is: 

g n (I)=I 0 e- A ' w (31) 

where the constant integration Io is the radiant intensity of the point source. Note that the 
above equation automatically satisfies the boundary condition: 

g m (~)=o 

[0069] By multiplying the RTE-SM by Q m for m=l . . .«>, the individual solutions can be 
15 superposed to get the final solution: 

I(T^) = |;c m g in (T)(L m . 1 (//) + L ra (//)). (32) 

m=l 

The values for coefficients Cm = 1 for all m, satisfy desired boundary conditions and they are 
shown to be accurate through real experiments as well as numerical Monte Carlo simulations. 
[0070] It is apparent that many advantages may be provided in the various embodiments of 

20 the systems and methods of the present invention. For example, in certain embodiments, 
weather, visibility, and range information corresponding to a light source captured in an 
image can be quickly determined. As another example, in certain embodiments, images can 
be altered to add or remove the appearance of multiple scattering effects. As yet another 
example, in certain embodiments, the techniques disclosed can be used for medical imaging 

25 of blood, tissues, or any other suitable portion of a body. As still another example, in certain 
embodiments, the techniques disclosed can be used for underwater imaging. As a further 
example, in certain embodiments, the results of the present invention can be used to narrow 
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the scope of the integral in a Monte-Carlo technique for modeling multiple scattering and 
thereby speed up that process. As a still further example, in certain embodiments, these 
techniques may be used for both absorbing and purely scattering media. As a yet further 
example, in certain embodiments, these techniques may be used to average the impact of 
multiple scattering over a long distance or a long time, and thereby minimize noise impacts. 
As a still further example, in certain embodiments, these techniques may be used as a first 
step before performing a more thorough Monte Carlo analysis of multiple scattering in the 
medium. 

[0071] Persons skilled in the art will thus appreciate that the present invention can be 
practiced by other than the described embodiments, which are presented for purposes of 
illustration and not of limitation, and that the present invention is limited only by the claims 
that follow. 
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